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We study the ground states of rotating atomic Bose-Einstein condensates with dipolar interactions. 
We present the results of numerical studies on a periodic geometry which show vortex lattice ground 
states of various symmetries: triangular and square vortex lattices, "stripe crystal" and "bubble 
crystal" . We present the phase diagram (for systems with a large number of vortices) as a function of 
the ratio of dipolar to contact interactions and of the chemical potential. We discuss the experimental 
requirements for observing transitions between vortex lattice groundstates via dipolar interactions. 
We finally investigate the stability of mean-field supersolid phases of a quasi-2D non-rotating Bose 
gas with dipolar interactions. 
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I. INTRODUCTION 

One of the most dramatic manifestations of the collec- 
tive quantum behaviour of a Bose-Einstein condensate 
is its unusual response to rotation. The formation of 
quantised vortex lines, around which the phase of the 
condensate wavefunction changes by 2n, is a direct con- 
sequence of the existence of a collective phase-coherent 
quantum state. The further ordering of quantised vor- 
tices at low temperatures into regular vortex lattices is a 
signature not just of Bose-Einstein condensation but of 
the presence of non-zero (repulsive) interactions between 
the constituent particles. 

Techniques for rotating condensates in ultra-cold 
atomic gases are now well established, and have led to 
very beautiful demonstrations of the physics of quan- 
tised vortices in these systems. Experimental studies of 
rapidly rotating atomic Bose gases have shown evidence 
for the formation of ordered arrays of very large num- 
bers of quantised vortices, P, 0, yj and detailed studies 
of the structure and dynamics of these vortex lattices 
have been performed. These studies are consistent with 
the expected properties for a Bose gas interacting with 
short range contact interactions (representing the van der 
Waals' forces) , for which the rotating groundstate is a tri- 
angular vortex lattice. 

The achievement of Bose-Einstein condensation of 
chromium-52Q opens up the possibility of experimen- 
tal studies of condensates with long-range interactions. 
Chromium has a large permanent magnetic dipole mo- 
ment, which leads to significant dipolar interactions in 
addition to the usual short-range interactions. It is of 
interest to ask how these long-range interactions affect 
the vortex lattice groundstates of the rotating conden- 
sate. It has been shownpj that the non-local interactions 
arising from dipolar interactions can lead to changes in 
the structure of the vortex lattices. The triangular vor- 
tex lattice can be replaced by vortex lattices of different 
symmetries: square lattice, "stripe crystal" and "bub- 



ble crystal" phases. (For a related study of square and 
stripe crystals see Ref0.) The vortex lattice structure 
is therefore a sensitive indication of the form of the in- 
terparticle interaction. However, the results of Ref. 
are restricted to the limit of weak interactions, when the 
mean interaction energy is smaller than the level spacing 
of the harmonic trap and the single particle states are 
restricted to the lowest Landau level (LLL) 0, @ . 

In this paper we extend the results of Ref . |3 beyond the 
weak interaction regime, and determine the vortex lattice 
structures to be found for arbitrarily strong interactions. 
For sufficiently weak interactions we recover the results 
of Ref. 13; we find groundstates that are the square lattice, 
"stripe crystal" and "bubble crystal" phases predicted in 
Ref. |3 in the LLL. We determine the effects of Landau 
level mixing on these states. We show that for strong 
interactions, these new vortex lattice states are unstable 
to collapse of the condensate. Our results establish the 
conditions required in order to observe transitions in the 
vortex lattice groundstate driven by dipolar interactions. 

We also describe the connections with possible "su- 
persolid" states of the non-rotating gas, and discuss the 
stability of these states. 

II. FORMULATION 

We consider a BEC of atoms with mass M which are 
confined in a harmonic trap with cylindrical symmetry 
about the z axis. We denote the trap frequencies by aiii 
and Wj_ in the axial and transverse directions, and the as- 
sociated trap lengths by a||j_ = y/WJ (MtJn±j. The atoms 
are taken to interact through both contact interactions 
and dipolar interactions, with a potential 

V(r) =gS 3 (r)+C d V d (r), (1) 

where g parameterizes the contact interactions, and Cd 
the strength of the dipole interaction potential, Vd(*°)- 
We consider the atomic dipole moments to be aligned 
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with the z axis, as in the recent experiments |4|, which 
sets the functional form for Vd(r). Since the classical 
dipole-dipole interaction can include also a delta-function 
contribution^], for clarity we define the dipolar interac- 
tion potential by 

V d (r^0)^ X l ±4 ~% (2) 
(x A + y l + z £ YI £ 

in which we exclude the point r = 0, such that any con- 
tact contribution of the dipole interaction is contained in 
the term proportional to g in Eq. We consider the 
ratio g/Cd to be a parameter that can be varied. For 
example in the chromium condensate g can be varied ex- 
perimentally by using a Feshbach resonance HE3 The 
system can remain stable even for a small attractive con- 
tact interaction (g < 0). (We note, however, that the 
q = Fourier transform of the effective potential, V 9= 0i 
should remain positive.) 

The results of Ref. 5 suggest that the mean field ap- 
proximation is excellent even for very dilute systems 
which may result from fast rotating gases. We thus em- 
ploy here the Gross-Pitaevskii model. We consider a sys- 
tem, rotating with angular frequency il, and containing 
a large number of vortices, such that near the center 
of the trap the particle density varies weakly over the 
scale of the vortex spacing. In the rotating frame, the 
Hamiltonian becomes equivalent to the Hamiltonian of a 
particle of c harg e q in a uniform magnetic field B, with 
qB/AI = 2fi[lll[!J. We discretise space and write the 
kinetic energy for this system in the form |13| 

£kin = -tJ2*t*je i4i3 '+cc, (3) 

where denotes the nearest neighbours of a square 

lattice with lattice constant a. In our calculations, 
the lattice constant a is chosen sufficiently small as to 
represent the continuum limit, with the particle mass 
M = ta 2 /h 2 . The "Peierls" phases are defined by the 



integral 

^ = ! ^ ' A-dl (4) 

calculated between the positions of the lattice sites i and 
j. We choose the Landau gauge for the vector potential 
A = B xy 1 which is convenient in our numerical calcula- 
tions. 

We determine the properties of a large array of vortices 
by studying the system in a periodic geometry. Impos- 
ing periodic boundary conditions under magnetic trans- 
lations requires that 

L x L y = 2ne 2 N v , (5) 

where iV v is the number of vortices in the simulation 
cell with dimensions L x ,L y . The magnetic length I = 
yJh/qB ~ y / h/2M£l arises as the natural length scale 
in this system. Since we consider a system contain- 
ing a large number of vortices, we have f2 ~ w i |l4l| so 
I ~ aj_/\/2. Vortex lattices in periodic systems have 
been studied previously for the case of contact interac- 
tions jiEj . 

Throughout this paper, we work in a quasi- two- 
dimensional (2D) regime, valid when the confinement in 
the longitudinal direction is sufficiently strong that the 
mean interaction energy (which is on the order of the 
chemical potential is small compared to the level 

spacing Tilh^ . We shall, however, allow mixing of Lan- 
dau levels for the in-plane motion, so shall work in the 
regime ~huj\\ ^> p, ~ hui± which implies a 2 /a\ <§; 1, i.e. 
that the trap is oblate. In the quasi-2D limit, p, <C hu>\\, 
the condensate wavefunction is a gaussian along the ro- 
tation axis, i.e., ^(x, y, z) — tp(x, y) e~ z ^ 2a n /(ir 1 ^^/ 2 ). 
Integrating the dipolar potential energy along the z-axis 
one obtains the effective dipolar interaction in two di- 
mensions 



V 2D (P) = I [V d (p, zi -Z 2 ) e 21 Z l g " dz 1 dz 2 ^-^^ e C 2 /4 [(2 + C 2 )Ko(C 2 /4) „ c2 ^ i(c 2 /4)]i (6) 

J J 7rajj 2v27rajj 
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where p = \fx* + y 2 , C = p/ a \\i an d K n are modified 
Bessel functions. When p is large compared to on, the 
expression © tends to l/p 3 , which is the dipolar energy 
for a 2D monolayer. We calculate l/p 3 by "Ewald sum- 
mation" for the present case of a periodic configuration of 
dipoles , and the remainder V^ 2D (p) — 1/ p 3 by direct 
summation within an area with radius several times the 
length aii (this combination is rapidly decaying with p so 
it converges quickly). 



For the energy minimization we use a variant of a 
norm-preserving relaxation algorithm which capitalizes 
on a virial relation |Is[ in order to fix the wavefunction 
norm. This is fed with an initial guess and is iterated 
until an energy minimum is reached. In order to investi- 
gate the vortex lattice groundstates, we initially simulate 
relatively large systems containing 7V V = 8 or 16 vortices. 
(Typical numerical grid sizes are N x ,N y ~ 50 — 60.) 
From these calculations, we find the range of values of 
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FIG. 1: The phase diagram for a trap asymmetry au/t = 0.4. 
The horizontal axis gives the ratio of the contact to dipo- 
lar interactions and the vertical axis is the chemical potential 
(normalized to the energy corresponding to the rotation fre- 
quency Q of the condensate). The condensate is unstable to 
collapse on the left of the solid line. Marked on the figure are 
the regions where vortex lattices of different symmetries are 
the ground state of the system. The lines are fits to numerical 
results marked by circles. 



the parameters over which the system is stable against 
collapse. We further identify the symmetries of the vor- 
tex lattice groundstates that appear, starting from differ- 
ent initial conditions to give unbiased searches of vortex 
lattices of different symmetry- types. We subsequently 
proceed to a systematic calculation of the energy for 
each of the above mentioned vortex lattice phases by 
using an appropriate simulation cell that is commen- 
surate with the translational symmetry of that phase. 
We use a square cell in order to obtain a square vor- 
tex lattice (e.g., N x = 10, N y = 10, N v = 1), a cell 
with N x /N y = V3 for a triangular vortex lattice (e.g., 
N x = 11, N y = 19, N v = 2), and accordingly for the 
other vortex lattice phases. The typical lattice spacing, 
implied by Eq. ©, is Ax = Ay = 0.251 We find the 
energy of each type of vortex lattice as a function of the 
parameter g/Cd, for various interaction strength values 
measured by the chemical potential \x. |l9j | 



III. VORTEX LATTICES 

The quantitative results depend on the asymmetry of 
the trap, aii/a± = a^/(\^2£). The phase diagram for 
the system au/£ = 0.4 is shown in Fig. ^ as a func- 
tion of g/Cd and of [i/hVL. The limit fj, <C h£l corre- 
sponds to the LLL. The region of stability of the system 
is shown by the solid line. Within the region of stability, 
we find a series of vortex lattice groundstates which in- 



clude the vortex lattices of the symmetries found in the 
LLL calculation 5]: triangular and square vortex lattices, 
stripe crystal and bubble crystal phases. (We refer the 
reader to Ref. |5j for illustrations of the particle distribu- 
tions in these phases.) The triangular vortex lattice is 
the lowest energy state for all values of g/Cd to the right 
of the dashed-dotted line. This includes, in particular, 
all positive values of g. On the left of the dashed-dotted 
line a square vortex lattice has energy lower than the 
triangular one. As g/Cd is further decreased the ground 
states are stripe crystal phases (in the region between the 
dotted and the dashed line) which can also be viewed as 
rectangular lattices of the vortices^. We solve for the 
ground states in this regime by varying the aspect ra- 
tio of the simulation cell and optimising the period of 
the stripes. As g/Cd is decreased the aspect ratio of the 
rectangular lattice varies, causing the separation between 
the stripes to become of larger period. For values of g/Cd 
to the left of the dashed line we find that a bubble crystal 
phase has the lowest energy. This is a triangular lattice of 
bubbles (clusters of particles), with the vortices arranged 
between the bubbles. The simplest bubble crystal phase 
has four double vortices per bubble which are arranged 
on the sites of the honeycomb lattice that is dual to the 
triangular lattice of bubbles 0. 

Our results, presented in Fig. ^ show that the values 
of g/Cd at which there are transitions between the vortex 
lattice phases are only weakly dependent on the chemi- 
cal potential, /i. Thus the transitions are closely given by 
the LLL theory; consistent with this the particle densities 
we find (not shown) closely resemble those in the LLLyj- 
We find that the main effect of going beyond this weak- 
interaction regime, is the appearance of the instability to 
collapse for large values of the chemical potential. The 
value of fi/hfl at which this instability sets in is a strong 
function of g/Cd- We associate this strong dependence to 
the varying particle density in these phases. The bubble 
crystal phase has a very inhomogeneous density, so, for 
a given average density, the maximum particle density is 
much larger in the bubble crystal phase than in the trian- 
gular vortex lattice. Indeed, in the LLL limit, fi <C hQ, 
the particle density at the center of the bubble diverges 
at the point of instability, which is g/Cd — > —8.49 for 
LLL calculation 5] applied to a±/£ = 0.4. [2(| 

Our results show that for large values of the chem- 
ical potential only the triangular vortex lattice is sta- 
ble. In this limit, fj, 3> hfl, the vortex cores are small 
compared to the vortex spacing, and the interactions be- 
tween the vortices arise from the kinetic energy of the 
flow around the vortices. This leads to a logarithmic re- 
pulsion between the vortices.|2l| and thus a triangular 
lattice groundstate. 

An estimate for the region of stability of the trian- 
gular vortex lattice in the strongly interacting regime, 
fi/HQ 3> 1, can be obtained by considering the condition 
for stability of the gas in the regions between the vortices. 
Denoting the spacing of the triangular vortex lattice by 
a v = [(47r) 1 / 2 /3 1 / 4 ]£, we look for an instability to modes 
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with wavevectors larger than q = a(2ir/a v ) of a homo- 
geneous Bose gas moving with local velocity (3h/(Ma v ) 
(a and j3 are numerical factors of order one). An analy- 
sis of the Bogoliubov spectrum (discussed in more detail 
below) shows that an instability occurs for 
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g+4vC d 

where A = \/3/(27r)[(27r) 2 a 2 - 2(3 2 ] and B = 
(27r) 2 3 1 / 4 a/v / 2. The functional form (jJJ) provides a good 
description of the boundary of stability of the triangular 
lattice in Fig. ^ at large (J./H&1 for reasonable choices of 
a and j3. This theory shows that, for large fi/hQ ^> 1, 
the boundary of stability of the triangular lattice tends 
to a fixed value g/Cd = — 47r + B-j- that increases with 
an/ 1. (To determine the precise value of this boundary 
would require a full solution of the Bogoliubov modes for 
the triangular vortex lattice in this regime.) 

We have investigated the effect of the trap asymmetry 
by repeating our numerical calculations for an ft = 0.8 
(a more symmetrical trap). The results for the region 
of stability, and the phase boundary between triangular 
and square vortex lattices are qualitatively similar to the 
case of Fig. ^ We find now that the transitions between 
the triangular and the square vortex lattices occur for 
larger values of g/Cd — —2.6 as compared to the case 
an 1 1 — 0.4. This is consistent with results in the LLL, 
and with the qualitative behaviour of J7J. Furthermore, 
the new vortex lattices are stable up to larger values of 
n/Ml (up to fi/Tifl ~ 12 for the square lattice) for this 
more symmetrical (less oblate) trap. 

In order to observe transitions in the vortex lattice 
groundstates induced by dipolar interactions, it is neces- 
sary both to tune g/Cd to a sufficiently negative value, 
using for example a Feshbach resonance, and also to en- 
sure that the particle density is sufficiently dilute (or 
interactions sufficiently weak) that /i/hQ is sufficiently 
small for the new vortex lattice states to be stable. We 
see in Fig. ^that this requires fi/Ml < 8 and g/Cd ^ 
—4.5 for the square lattice to be observed. (Stability 
of the stripe crystal and bubble crystal requires more 
stringent conditions.) Our numerical calculation gives 
directly the number of atoms per vortex that correspond 
to the above conditions. For the parameters relevant 
for 52 Cri, and choosing Q — 2n x 100Hz, the above 
conditions require that the number of particles per vor- 
tex be v = n 2 d/n v < 1300. Similarly, for a\\/£ = 0.8, 
the conditions for stability of the square vortex lattice 
(fi/nn < 12, g/Cd < -2.5) lead to v < 3000. While it 
is challenging to achieve a vortex lattice in a sufficiently 
dilute atomic gas, we note that these conditions are in 
well excess of what has been achieved in rapidly rotating 
Rubidium condensates. 0, Q 

Previous calculations of the effects of dipolar interac- 
tions on vortex lattices performed in a trap geometrv|22j 
found direct transitions from triangular vortex lattice to 
collapse. It is not straightforward to make direct connec- 
tions between the present results for the phase diagram 



of a homogeneous vortex lattice, with those studies on in- 
homogeneous trap ped systems. However, we note that in 
the studies of Ref . |23 the number of vortices is relatively 
small, so the confinement induces a sizeable variation in 
the particle density even on the scale of a vortex spacing. 
It is possible that, owing to the higher particle density 
at the centre of the trap, a trapped system with a small 
number of vortices is less stable to collapse. Our results, 
Fig-d apply to a system with a large number of vortices 
such that density variation over the scale of the vortex 
lattice constant is small. 



IV. NON-ROTATING DIPOLAR BOSE-GAS 

We now turn to investigate the non- rotating gas. It 
is natural to ask whether the stripe and bubble states 
that appear for a rotating condensate are manifestations 
of a "supersolid" , in which density-wave order appears 
spontaneously, in the non-rotating groundstate. In the 
latter case one would expect its response to rotation to 
involve the vortices being placed in regions of low particle 
density. Depending on the symmetry of the supersolid 
state this could lead to particle distributions similar to 
those of the stripe crystal and bubble crystal states. 

A further motivation in this direction if offered by the 
development of a roton minimum in the excitation spec- 
trum of a homogeneous condensate with dipolar inter- 
actions [2^, |3 which can touch zero energy at a finite 
wavevector. This indicates an instability of the homo- 
geneous Bose gas and is suggestive of a possible phase 
transition to a density- wave state (e.g. a condensate of 
Bogoliubov modes with non-zero wavevector). 

The Bogoliubov modes of a quasi-2D homogeneous 
condensate with dipolar interactions were discussed in 
Ref. |2^. For a gas with density n 2 d the chemical po- 
tential is fi — ri2dVq=o, and the excitation spectrum is 

E (l) = \ e q( e q+ n 2dV q J , where e q = h 2 q 2 / (2M), and 



V q = 



47rC d 



2na 



2TrC d qe q a n /2 erfc(ga 



(8) 



is the Fourier transform of the effective 2D interaction 

m. 

If, as above, one considers g to be tuned to negative 
values, the Bogoliubov mode develops a roton minimum 
that becomes unstable eve n for the quasi-2D regime, fj, -C 
fiu}\\, for g/Cd = — 47r + V2ir 3 Cdn -2d / ( a\\Tiu\\) ■ (The case 
g < was not considered in Ref. I25L and consequently 
the instability discussed there, for g > 0, occurs outside 
of the regime of validity of the quasi-2D approximation, 
in the crossover to the 3D regime /i ~ ?La;||.) 

We have looked numerically for supersolid states in 
the quasi-2D Bose gas. We set f2 = in the code used 
for the numerical calculations described above, chose pa- 
rameters close to the instability, and looked for stable 
(or metastable) states of the system. For all parame- 
ter values that we have tested (which are in the quasi-2D 
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regime) , we find that the only stable state is the homoge- 
neous Bose gas. Configurations with "supersolid" density 
wave order (showing a triangular lattice of peaks in the 
density) are unstable in all cases that we have studied. 
The collapse mechanism involves a process in which all 
particles accumulate into one of the peaks, generating a 
region of the sample in which the density diverges. We 
have tested this conclusion by independent numerical cal- 
culations (based on a discretisation in momentum space, 
and minimisation by the conjugate gradient method) and 
have found consistent results. Thus we conclude that, 
within mean-field theory, there are no stable supersolid 
phases for an infinite quasi-2D Bose gas with interaction 

V. CONCLUSIONS 

We have determined the phase diagram for a rapidly 
rotating atomic Bose gas with dipolar interactions, be- 
yond the weak interaction regime. Our results estab- 
lish that dipole interactions can drive transitions from 
the triangular vortex lattice into vortex lattices of vari- 
ous different symmetries. These provide estimates of the 
experimental parameters that are required in order to 
observe vortex lattice groundstates of types other than 
triangular. 

We also investigated the fate of a non-rotating quasi- 
2D dipolar Bose gas at the point at which the homoge- 
neous Bose gas becomes unstable to Bogoliubov modes 



with a finite wavelength. Our calculations indicate that, 
within the mean-field approach used, there do not exist 
stable supersolid states, but that the homogeneous Bose 
gas becomes unstable to collapse. 

The stripe and bubble phases discussed above for a 
rotating gas are therefore not manifestations of stable 
supersolid phases of the non-rotating gas. Rather, one 
should view the rotation as acting to stabilise the density- 
wave order that the dipolar interactions tend to induce in 
the dipolar Bose gas. The rotation introduces a length- 
scale, the magnetic length £, which, in the weakly inter- 
acting regime /i/htt -C 1, sets the minimum wavelength 
over which the condensate wavefunction can vary This 
minimum lcngthscale prevents collapse of the system pro- 
vided interactions are sufficiently weak. 

Note added: While preparing this work for publica- 
tion we learned that G. Shlyapnikov and P. Pedri have 
reached the same conclusion: that putative supersolid 
states in the quasi-2D Bose gas with the interaction 
are unstable. 1271 
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